SUBROUTINE insertPints(npt,maxpts,ymax,y,z)
IMPLICIT NONE
INTEGER, PARAMETER :: prec=KIND(1.0D0)    ! defines precision of problem, # of decimal points
INTEGER, PARAMETER :: iprec=SELECTED_INT_KIND(9)    ! defines precision of integers
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! purpose: gets boundary condition, stin, at time tt
!! parameters:
INTEGER(iprec),INTENT(INOUT) :: npt
INTEGER(iprec),INTENT(IN) :: maxpts  ! sediment solution type
REAL(prec),INTENT(IN) :: ymax  ! max horizontal distance between points
REAL(prec),INTENT(INOUT),DIMENSION(maxpts) :: y,z  ! cross section points
!! error conditions:    returns iflag /= 0
!! side effects:        none known
!! bugs:                none known
!! author:              Jianchun Victor Huang 12-17-2009, USBR
!! references:
!! notes:
!---------------------- local declarations ------------------------------------
INTEGER(iprec) :: k,nInsert,kk,nptNew
REAL(prec) :: y1,y2,z1,z2,dy,dz
k=1
DO while (k<npt)
    y1=y(k); y2=y(k+1)
    IF (abs(y2-y1)>yMax) THEN
        ! insert points
        nInsert=ABS(y2-y1)/dmax1(yMax,0.1)
        nptNew=npt+nInsert
        IF (nptNew>maxpts) THEN
            write(*,*) "ERROR in insertPints: too much points"
            write(*,*) "maxpts=",maxpts
            write(*,*) "npt=",nptNew
        ENDIF
        ! shift nInsert points
        DO kk=nptNew,k+1+nInsert,-1
            y(kk)=y(kk-nInsert)
            z(kk)=z(kk-nInsert)
        ENDDO
        ! interpolate nInsert points
        z1=z(k); z2=z(k+1)
        dy=(y2-y1)/real(nInsert+1)
        dz=(z2-z1)/real(nInsert+1)
        DO kk=1,nInsert
            y(k+kk)=y1+real(kk)*dy
            z(k+kk)=z1+real(kk)*dz
        ENDDO
        npt=nptNew
        k=k+nInsert+1
    ELSE
        ! skip to next point
        k=k+1
    ENDIF
ENDDO

END SUBROUTINE insertPints


SUBROUTINE insertPointAtY(npt,maxpts,yp,y,z)
IMPLICIT NONE
INTEGER, PARAMETER :: prec=KIND(1.0D0)    ! defines precision of problem, # of decimal points
INTEGER, PARAMETER :: iprec=SELECTED_INT_KIND(9)    ! defines precision of integers
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! purpose: add a point at yp
!! parameters:
INTEGER(iprec),INTENT(INOUT) :: npt
INTEGER(iprec),INTENT(IN) :: maxpts  ! sediment solution type
REAL(prec),INTENT(INOUT) :: yp  ! y coordinate with poinit added
REAL(prec),INTENT(INOUT),DIMENSION(maxpts) :: y,z  ! cross section points
!! error conditions:    returns iflag /= 0
!! side effects:        none known
!! bugs:                none known
!! author:              Jianchun Victor Huang 10-25-2010, USBR
!! references:
!! notes:
!---------------------- local declarations ------------------------------------
INTEGER(iprec) :: k,nInsert,kk,nptNew
REAL(prec) :: y1,y2,z1,z2,dy,dz
pt: DO k=1,npt-1
    y1=y(k); y2=y(k+1)
    IF(abs(y1-yp)<=0.1) Then
        y(k)=y1
        nptNew=npt
        EXIT pt
    ELSEIF(abs(y2-yp)<=0.1) Then
        y(k+1)=y2
        nptNew=npt
        EXIT pt
    ELSEIF(y1<yp .and. y2>yp) Then
        nptNew=npt+1
        DO kk=nptNew,k+2,-1
            y(kk)=y(kk-1)
            z(kk)=z(kk-1)
        ENDDO
        z1=z(k); z2=z(k+1)
        y(k+1)=yp
        z(k+1)=z1+(yp-y1)*(z2-z1)/(y2-y1)
        exit pt
    ELSE
        nptNew=npt
        cycle pt
    ENDIF
ENDDO pt
npt=nptNew
END SUBROUTINE insertPointAtY